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We present a new numerical technique for 
modeling the flow around multiple ob- 
jects moving in a fluid. The method tracks 
the dynamic interaction between each 
particle and the fluid. The movements of 
the fluid and the object are directly cou- 
pled. A background mesh is designed to fit 
the geometry of the overall domain. The 
mesh is designed independently of the pres- 
ence of the particles except in terms of 
how fine it must be to track particles of a 
given size. Each particle is represented 
by a geometric figure that describes its 
boundary. This figure overlies the mesh. 
Nodes are added to the mesh where the 
particle boundaries intersect the back- 
ground mesh, increasing the number of 
nodes contained in each element whose 
boundary is intersected. These additional 
nodes are then used to describe and track 
the particle in the numerical scheme. Ap- 
propriate element shape functions are de- 
fined to approximate the solution on the el- 
ements with extra nodes. The particles 



are moved through the mesh by moving 
only the overlying nodes defining the 
particles. The regular finite element grid 
remains unchanged. In this method, the 
mesh does not distort as the particles move. 
Instead, only the placement of particle- 
defining nodes changes as the particles 
move. Element shape functions are up- 
dated as the nodes move through the ele- 
ments. This method is especially suited 
for models of moderate numbers of moder- 
ate-size particles, where the details of 
the fluid-particle coupling are important. 
Both the complications of creating finite 
element meshes around appreciable num- 
bers of particles, and extensive remesh- 
ing upon movement of the particles are 
simplified in this method. 
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Glossary 



^ 


gravitational acceleration 


t/co 


in. 


mass of particle 


x,y 


n 


unit normal vector 


e 


P 


pressure 


n 


r 


radius of cylinder 


M 


R 


radial distance from center of cylinder 


<P 


T 


temperature 


p 


u 


velocity vector 


cr 


Mx 


velocity in x direction 


V 


Uy 


velocity in y direction 


V- 


u 


magnitude of velocity 


V^ 



constant velocity far from particle 

rectangular Cartesian coordinates 

penalty parameter 

angular velocity 

viscosity 

stream function 

density 

shear stress 

gradient operator 

divergence operator 

Laplacian operator 



77 



Volume 103, Number 1, January-February 1998 

Journal of Research of the National Institute of Standards and Technology 



1. Introduction 
1.1 Background 

Fluid-solid systems are classified based on the size 
and number of particles (or other objects). Particle sizes 
range from aerosols to large objects. The focus of an 
analysis may be the detailed movement of a single par- 
ticle, or the interaction of a large number of objects. 
Different numerical techniques are appropriate for dif- 
ferent classes of problems. In systems with very small 
particles, such as aerosols, the number of particles is 
typically very large. In these systems the effect of the 
particles on the flow is often ignored. Corrections can 
be made to account for the momentum transferred to the 
particles, but there are too many particles to consider 
the details of the coupled dynamical interaction be- 
tween the fluid and the particles. 

Details of the fluid-particle interactions are important 
in systems such as fluidized beds, steam turbines, natu- 
ral gas pipelines, and sprayers. In such systems it is 
important to know how particles interact to affect their 
trajectories, which is seen most clearly when particles 
make contact with one another. A simulation must cou- 
ple the dynamics of the fluid and the particles to achieve 
an accurate solution. 

One approach for such systems is to mesh the fluid 
surrounding the particles, taking into account the inte- 
rior boundaries created by the particles. Such meshes 
are necessarily complex because of the large curvature 
of the particle boundary compared with the boundaries 
of the larger domain. Consequently, it is convenient to 
use unstructured meshes. The ability to use unstructured 
grids is a potential advantage of finite element methods 
compared to finite difference techniques. When the par- 
ticle or particles move, the mesh is allowed to distort to 
some degree. This distortion tends to reduce the accu- 
racy of the simulation. After the particle(s) travel some 
distance the distortion becomes too great, and either the 
entire domain or a region around each particle is 
remeshed. 

Tezduyar and collaborators [1-5] are developing ef- 
ficient algorithms for moving the particles and updating 
the mesh. They have developed what they call a De- 
forming-Spatial-Domain/Space-Time procedure for 
simulating moving particles in two or three dimensions. 
Triangular (or tetrahedral in 3D) elements are used to 
discretize the fluid phase surrounding the particles. Spe- 
cial narrow elements are used near each particle to track 
the boundary layer. The elements near each particle 
deform as the particle moves through the mesh. 
Remeshing is necessary when the elements become 
overly distorted. Both remeshing and checking the 
mesh for distortion can be computationally intensive 
processes if they are done at each time step. Tezduyar 



et al. discuss that this is particularly true for computa- 
tions on parallel-architecture machines because of the 
time required to apportion the meshing tasks among the 
processors. Consequently, remeshing is typically per- 
formed only at specified timesteps. This creates a po- 
tential for increased error if too much distortion occurs 
between remeshing. As a result, it appears that this 
technique is most suitable for slowly deforming spatial 
domains. One means they employ to limit the extent to 
which remeshing is necessary is to translate the entire 
mesh at the velocity of the center of mass of the set of 
spheres. 

Hu, Joseph, and Crochet [6] describe a method for 
direct simulation of flowing "particles" in two dimen- 
sions, using unstructured triangular finite element 
meshes. At each time step they remesh, generating a 
new mesh that bears no straightforward relationship to 
the previous mesh. This then requires the regeneration 
of element-connectivity information, which is neces- 
sary for the finite element method to know how the 
solution interacts between adjacent elements. This is 
facilitated by a new mesh-data structure that can be 
processed more efficiently to find the position of a 
given node. The solution from the previous time step 
must be interpolated onto the new mesh. Their method 
can track computational domains that move drastically 
with time. They use a stable time-integration scheme in 
which "at each time step, the positions of the particles 
are updated explicitly, the computational domain is 
remeshed, the solution at the previous time is mapped 
onto the new mesh, and finally the nonlinear Navier- 
Stokes equation and the implicitly discretized Newton's 
equations for particle velocities are solved on the new 
mesh iteratively" [6]. The difficulties they address are 
remeshing at each time step, the interpolation of the 
solution onto the new mesh, and finally the general 
issue of solving the system in a stable and efficient 
manner. 

A significant issue of both of these techniques is the 
establishment of robust criteria for acceptable vs unac- 
ceptable distortion between remeshing. Remeshing can 
be computationally intensive, particularly in three di- 
mensions, and especially for parallel computations. Er- 
ror is introduced every time the solution at a given time 
step must be projected onto a new mesh. While ap- 
proaches involving unstructured meshes with remesh- 
ing have proved to work well for at least certain sets of 
problems, we were motivated to seek a technique that 
would avoid the computational difficulties of remesh- 
ing, and interpolation of the solution onto new meshes. 

Our method was motivated by that of Fogelson and 
Peskin [7]. They developed a fast finite difference 
method for solving Stokes' equations in three dimen- 
sions to describe the motions of elastically deformable 



78 



Volume 103, Number 1, January-February 1998 

Journal of Research of the National Institute of Standards and Technology 



suspended particles (rigid particles are modeled as 
slightly deformable). The domain is discretized on a 
regular lattice. They represent each particle by a sepa- 
rate set of points overlying the lattice. Equations of fluid 
motion are applied at all points on the lattice, including 
those overlain by particles. The particle nodes move at 
the local fluid velocity. However, the local fluid velocity 
is affected by the presence of the particle through a 
force-density term that accounts for the resistance of the 
particles to deformation. It affects the fluid near the 
particle. The particles are modeled as regions of the 
fluid with elastic cohesive forces. This allows them to 
solve the entire domain as a single phase, with no fluid- 
solid boundaries introduced by the particles. The result 
is a very fast technique at the expense of some accuracy 
about the fluid-particle interaction. In addition to the 
internal elastic links, any external forces are also han- 
dled in a straightforward way, including interparticle 
forces. Since this is a purely Stokes model, it does not 
account for the inertia of the particle or the fluid. 

Unverdi and Tryggvason [8] briefly survey different 
methods of modeling sharp fronts that illustrates the 
breadth of different techniques that are available. In this 
paper they use an approach that is somewhat similar to 
that of Fogelson and Peskin. However, it tracks the inter- 
face explicitly, which allows them to consider different 
phases as distinct. To look at the motion of bubbles in a 
fluid, they use a regular finite difference grid onto 
which they overlay an unstructured grid that represents 
the particle surface. The interface region between the 
particle surface and the surrounding fluid is given a 
thickness, on the order of the mesh spacing, to allow a 
more gradual transition between the phases. This is done 
to provide stability and smoothness. The main advan- 
tage of this method is its ability to track multiple inter- 
facial boundaries in the same region of the grid. 

1.2 A New Technique 

Our work follows a somewhat different approach. We 
wanted to avoid the difficulty of remeshing a significant 
fraction of the domain. This led us to an approach using 
overlying sets of nodes to represent the particles. We 
also wanted to be able to render the solution around a 
particle as accurately as possible. In adapting this ap- 
proach to finite elements it appeared natural to make the 
extra nodes belong to the elements in which they appear. 
We have developed a dynamic finite element type in 
which the shape functions used on a given element are 
allowed to change in time as a particle boundary moves 
into and out of an element. 

A background mesh is designed to fit the geometry of 
the overall domain. The mesh is designed independently 
of the presence of the particles except in terms of how 



fine the mesh must be to track particles of a given size. 
The particles are represented by circles which overlie 
the mesh. Nodes are added to the mesh where the parti- 
cle boundaries intersect the (background) mesh, increas- 
ing the number of nodes contained in each element 
whose boundary is intersected. These additional nodes 
are then used to describe and track the particle in the 
numerical scheme. Appropriate element shape func- 
tions are defined to approximate the solution on the 
elements with extra nodes. The particles are moved 
through the mesh by moving only the extra, overlying 
nodes defining the particles. Otherwise, the mesh re- 
mains unchanged: the background mesh does not distort 
as the particles move. Element shape functions are de- 
pendent on the placement of the particle nodes with 
respect to the element sides, and change dynamically as 
the particles move through the mesh. Shape functions 
are updated as the particle nodes move through the 
elements. We demonstrate our technique with flow ex- 
amples for which exact solutions are available. 

A primary advantage of this method is that there is no 
mesh distortion as particles move. This is significant 
because remeshing tends to be time-consuming and er- 
ror-prone. We also avoid extensive checking for whether 
remeshing is necessary. The particles intersect underly- 
ing elements at different locations as they move, and the 
element shape functions change to account for this 
movement. The only additional storage required is the 
current location of the particle nodes, and their relative 
position along an element side. At each time step of a 
transient analysis, the translation and rotation of each 
particle are found directly from the continuum equations 
for the fluid, and Newton's equations for the particles, 
and the nodes that represent that particle are moved 
accordingly. Specifically, intersection points where par- 
ticle boundaries meet element edges, and variables de- 
scribing positions along element edges are updated. The 
method is easily extended to particles with irregular 
shapes, and different shapes for different particles. 

The various methods described above are best suited 
to different classes of problems. The first approach, in 
which the mesh is allowed to distort, is well suited to the 
accurate solution of a small number of particles. With 
larger numbers of particles there is a temptation to 
remesh less often because of the expense, which com- 
promises accuracy. The finite difference approach using 
an overlying mesh is well suited for large numbers of 
particles, especially when the details of the flow near 
the particles are not as important. Our method avoids 
the complications of creating finite element meshes 
around appreciable numbers of particles, and the exten- 
sive remeshing required to track the movement of parti- 
cles is eliminated. At the same time, it incorporates the 
particle boundaries directly into the mesh, allowing for 
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an accurate solution near the particles. This tends to 
make this technique more appealing for a system with a 
moderate number of particles, which cover most of the 
domain. 

If large areas of the domain contain no particles, then 
computational resources are unnecessarily concentrated 
there. Nonetheless, we believe that this technique, or an 
extension of it, may also prove to be useful for systems 
involving smaller numbers of particles. A possible ex- 
tension is simply to add extra "rings" or layers of nodes 
around the moving objects as needed to accurately rep- 
resent the solution, so that the node density is increased 
primarily only in localized regions around the particles. 
This would also be very useful for modeling flows at 
higher Reynolds numbers, where the local fluid-particle 
interactions are more complex. Another possibility is to 
refine the mesh locally around each particle in a pre- 
scribed way. These well-defined "refinement zones" 
would move with the particles, but the background mesh 
would remain undistorted. 



2. Equations of Motion 
2.1 Flow Equations 

The equations describing conservation of momentum 
and mass are as follows. (For definitions of the various 
quantities, see the glossary.) 

Momentum equation: 

Po(-^ + M ■ Vh) = - V/J + ^tV^H + pg, (1) 

where u = (u^, Uy, mJ. 
Continuity equation: 

V ■ M = 0. (2) 

A standard finite element methodology for incompress- 
ible flow is to replace the continuity equation with a 
"penalty" constraint (as described in Hughes et al. [9]), 
in order to eliminate the need to solve for the "pres- 
sure" field directly. This is often useful, as long as the 
particle boundary conditions are not direct functions of 
pressure. 

Penalty equation: 



V ■ M = 



ep. 



(3) 



Our code also has other capabilities, such as modeling 
of heat transfer and chemical reactions, that are of inter- 
est in a fluid-particle simulation. These capabilities are 
easily extended to the particles and particle boundaries. 
However, they are not germane to the present discussion. 



so for clarity we omit them. Boundary conditions for the 
particles prescribe no slip at the fluid-particle interface. 
Boundary conditions: 



Mfiuid = Msoiid at particle boundaries. 



(4) 



Presently, we are considering only solid particles, but 
this technique should be readily extensible to fluid drops 
or bubbles, yielding a two-phase fluid-fluid system. 

2.2 Equations to Represent Particle Movement 

For the purposes of developing and testing this tech- 
nique we are presently working only on two-dimen- 
sional flow. This is a constraint of practicality for devel- 
opment. It is straightforward to extend the method to 
three-dimensional flows. The movement of a 2-D "par- 
ticle' ' is determined by three degrees of freedom: two 
translational and one rotational. All nodes defining an 
individual particle are constrained to have the same 
translational and rotational movement. Components for 
rotation vary according to their position on the circle 
representing the boundary of the particle: each has the 
same speed, but a different direction of movement in a 
Cartesian sense. The net velocity of each particle node 
is a sum of these contributions. An additional constraint, 
conservation of angular momentum, is added to the set 
of global equations for each particle corresponding to 
the rotational degree of freedom. Particles conserve 
both linear and angular momentum. 

Conservation of linear momentum: 

rripg + I (n ■ (t)As = rripi ^'^ ^ ). (5) 



Conservation of angular momentum: 



/ 



(r X (« ■ a)) = (^)I. 



(6) 



To evaluate these constraints on each particle, we define 
one-dimensional line elements along the particle 
boundary. Each line element lies within a two-dimen- 
sional four-node element of the background mesh. 
Within this element is constructed a separate two-di- 
mensional four-node "sub-element" associated with 
the particular line element. This new sub-element is 
handled exactly like any other four-node element, but is 
formed solely for the purpose of computing normal 
gradients of the velocity along the line element, for 
computation of the stress boundary condition. Two of 
the corner nodes of this sub-element are the end-points 
of the line element, which are the intersections of the 
particle boundary with the particular background ele- 
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ment. The other two corner nodes are the adjacent cor- 
ner nodes of the background element. Thus, two of the 
sides of the sub-element lie along sides of the back- 
ground element, one side approximates the boundary of 
the particle that overlaps the background element, and 
the fourth side of the sub-element is a diagonal of the 
background element. This construction is illustrated in 
Fig. 1. 




Fig. 1. Construction of a sub-element for the computation of a normal 
to a particle. 



3. Variable Finite Element Shape 
Functions 

As discussed above, we are currently limited to two- 
dimensional simulations. We are using bilinear 
(straight-sided) 4-node quadrilateral elements for the 
background mesh. We superpose on this ' 'background' ' 
mesh circles that represent the outlines of particles. This 
is generalizable to somewhat irregular shapes. Repre- 
sentation of highly contorted shapes will require further 
development. We add an extra node at each intersection 
of a circle with an element side. For a given 4-node 
element, we consider cases where the particle intersects 
the elements edges at either 2 or 4 points. It is very rare 
for the particle to intersect an element at exactly 1 or 3 
points, because the time dimension is finite differenced, 
and consequently, discrete. Intersection at 1 or 3 points 
would require the boundary of the particle to be exactly 
tangent to an element side. Numerically the particles 
move in small increments, so the chances of the edges 
coinciding tangentially to within machine accuracy are 



minute. Currently, we simply ignore such points, assum- 
ing they do not coincide. As yet, we have not observed 
exact coincidence to occur. 

We developed shape functions for the resultant 6-node 
or 8-node elements, based on the method described in 
Zienkiewicz [10] for generating "serendipity" ele- 
ments. In the finite element method, a solution to the 
transport equations for a given dependent variable is 
represented over an element by summing the contribu- 
tion from each shape function, corresponding to each 
node, over all of the nodes of the element. Lagrangian 
interpolation functions are generally used as the shape 
functions because they yield the convenient property 
that the nodal values of the dependent variables are 
available directly, without summing contributions. For 
example, temperature on an element, T, would be repre- 
sented by: 



T=^(i>,Ti. 



(7) 



Ti factors are coefficients (constant at a given timestep) 
corresponding to the temperature at each node i of the 
element, and the <^, are the nodal basis functions. 

Following standard finite element methodology, each 
element on the "physical" domain (the mathematical 
representation of the physical domain), on which the 
problem is posed, is mapped locally onto a regular 2X2 
square element, called the "parent" element, centered 
around the point (0,0) (see Fig. 2). This allows the use 



(■1>-1) 



(1,1) 



Fig. 2. The parent element, onto which all elements in the physical 
domain are mapped. 
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of one set of standard shape functions, which is defined 
on this unchanging ' 'parent' ' . The representation of the 
solution on each element is actually formulated on this 
parent element. We designate a Cartesian coordinate 
system on this parent element as {s, t]. Both s and t 
range from — 1 to 1 . In this system, the shape functions 
for a bilinear 4-node element, where the nodes are num- 
bered counterclockwise starting at the lower left, are 
defined below. 

Standard bilinear finite element shape functions for a 
quadrilateral element are: 



iV,=^(i-^)(i-0, 



N2 = ^{\+s){l-t), 



(8) 



(9) 



mapped domain and the "physical" element on the 
unmapped domain. It also eliminates the need to sum 
the polynomials representing the solution on each ele- 
ment. The computation of the nodal solutions is direct 
(see Zienkiewics [10].) 

For the new, 6-node element, two new shape func- 
tions are needed corresponding to the two new nodes. 
Consider, for example, the case with a particle node 
along the bottom (side 1 ) and another on the right side 
(side 2). These new nodes are designated nodes 5 and 6, 
respectively. Letp represent the s value of the extra node 
along side 1 , and q represent the t value of the extra node 
along side 2. The new nodes are located a.t(p, — 1) and 
(1, ^). The shape functions for each of the extra nodes 
is piecewise continuous. 

New bilinear finite element shape functions for parti- 
cle nodes are: 



N, = -{l+s){l+t), 



N, = -(l-s)(\+t), 



(10) 
(11) 



Each function has a value of 1 at the node for which it 
is defined, and a value of at all other node points. This 
is the property that allows the solution at the nodes to be 
obtained directly. When new nodes are added, care must 
be taken to adjust the shape functions at all affected 
nodes to account for this. Otherwise, unintended contri- 
butions are erroneously added to the representation, dis- 
torting the solution. Again, this is standard finite ele- 
ment methodology. The novelty lies in the use of the new 
elements to track the particle boundaries. 

For elements that contain extra nodes, for each addi- 
tional node an additional bilinear shape function is 
added to the representation of the solution on that ele- 
ment. The specific shape function that is added for each 
extra node is based on the position of that node along the 
side, and depends on which side of the element it is on. 
Shape functions for the end nodes (corner nodes) of 
sides containing an extra node must be modified. Addi- 
tional terms are added to them to ensure that their shape 
function value is zero at the new nodes as well as the 
other "original" (here, corner) nodes in the element. 
This is an important facet of the finite element approach 
based on Lagrangian polynomials. The requirement is 
that at each node of an element only one shape function 
contributes, the one corresponding to that node, and the 
value of the shape function is 1 . With this formulation, 
the value of a dependent variable at a node is simply the 
value of the coefficient of the shape function for that 
node. This avoids any need to compute an inverse map- 
ping between the "parent" element on the locally 



Ns = < 



Ne 



(i-m + s) 

2(1 +p) 

(1 - 0(1 - s) 
2(1 -p) 

(l+s)(l+t) 
2(1 +q) 

(l+s)(l-t) 
2(1 - q) 



fOTS^p 



for s > p. 



foxt^q 



foT t> q. 



(12) 



(13) 



At points 5 and 6 the values of shape functions A^i — A'4 
must equal 0. Consequently, these four shape functions 
must be changed according to the location of the extra 
node. 

The modified bilinear finite element shape functions 
for corner nodes of a particle element are: 



N{ = N:-^Ns(l-p), (14) 



m = N2-^N,(l+p)-^N6(l-q), (15) 



m = N,--N6(l+q), (16) 



N4 = N4 (unchanged) . (17) 

The entire set of six shape functions on the new 
parent element is: 
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N>=l{l-s){l-t)-^m{l-p), (18) 



N2 = \{1+ s){l - t) - I A^5(l +P)- ^N^H - q). 



(19) 



N, = ^(l+s)(l+t)-^Ne(l+q), (20) 



N4 = ^(l -^)(l+0, 



(21) 



N5=< 



N6=< 



{l-t){l+s) 
2(1 +p) 

(1 - 0(1 - s) 
2(1 -p) 

(l+^)(l+0 
2(1+?) 

(l+^)(l-0 

2(1 - q) 



fOTS^p 



fOT S> p, 



for t: 



for t> q. 



(22) 



(23) 



We use a subparametric mapping from each ' 'particle 
element' ' (each element containing part of a particle) in 
the physical domain to the parent element: more nodes 
are used to approximate the solution on a particle ele- 
ment than are used to define the mapping. This is done 
to preserve a linear mapping, to avoid unnecessary dis- 
tortion in moving from the physical domain to the lo- 
cally mapped domain. This is sensible, since in this 
technique the movement of the particles does not cause 
any distortion of the elements of the background mesh. 
Only the four corner nodes are necessary to define the 
Jacobian of the transformation from the coordinate sys- 
tem on the element in physical domain to (s, t) coordi- 
nates on the parent. The values of p and q are updated 
following each step in a transient analysis with moving 
particles. There are significantly fewer steps involved in 
updating the mesh than in remeshing schemes. (See, for 
example, Johnson et al. [1] and Hu et al. [6] for discus- 
sions of particular mesh-updating strategies.) 

For cases in which an edge of a particle moves very 
close to an element side, there can be two intersections 
of the particle and the element along the same edge. In 
that case a sUghtly different set of linear shape functions 
is defined for that element. Suppose there is a resulting 
6-node element with 2 extra nodes along side 1 whose 



s coordinates have values of p and q. Assuming/? is less 
than q, let pi = (p — l)/2, qi = (q + l)/2. The corre- 
sponding six shape functions are: 



M = ^ (1 - ^)(1 - t) -^N,(l -p)- ^Ne(l - q). 



(24) 



N2 = \(l+s)(l-t)-^N,(l+p)-^N,(l+q), 



N, = -(l+s)(l+t), 



N, = ^(l-s)(l+t), 



(25) 
(26) 

(27) 



A^5=< 



(s + l)(t - 


1) 


r- 2(p+l) 




(S - q)(t - 


1) 


2(p-q) 




(S - q)(t - 


1) 


A(q - 1) 




^ (1 - S)(t - 


1) 



A(q - 1) 



A^6=< 



(s + l)(t - 1) 


r- 4(p+l) 


(s-p)(t-l) 


4(p+l) 


(s-p)(t-l) 


2(p-q) 


K, (s- m - 1) 



2(q - 1) 



fOT S <p 

foT p <s <q 
for q< s < q\ 
for q\<s, 

for s <p\ 
foT Pi <s <p 
foT p <s <q 
foT q <s. 



(28) 



(29) 



(A^5 and N^ correspond to the two extra nodes on side 1.) 
We have incorporated this moving-boundary method 
into our own software, which we wrote in C++ for 
simulating and analyzing fluid flow with associated 
heat and mass transfer. It was written with an object- 
oriented design in order to facilitate adding new fea- 
tures, such as new moving-boundary methods, new sets 
of global equations, new boundary conditions, etc. For 
additional information about our software, see Peskin 
and Hardin [11]. 
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All of the simulations presented in this paper are 
two-dimensional representations of flows. We are cur- 
rently working on a three-dimensional version of this 
method, which is a straightforward extension of the two- 
dimensional case. Eight-node brick elements are used 
for the finite element mesh. A single element can be 
totally enveloped in a spherical particle, or partially 
enveloped. When only a portion of a brick element con- 
tains part of a sphere, the sphere intersects that element 
along either 3 or 4 of the brick edges in the most com- 
mon cases. Analogous to the 2-D cases, there are other, 
less-common possibilities that should also be accounted 
for. Different sets of shape functions are defined for 
elements corresponding to different types of sphere-par- 
ticle overlap. 



4. Test Problems 

We present simulations to test this moving-boundary 
technique. We first compare our numerical solutions to 
exact solutions for a set of Stokes-flow problems. We 
then present a solution of Navier-Stokes flow around 
multiple particles in a channel. 

Flow patterns around circular two-dimensional ' 'par- 
ticles' ' represent cross-sections of flows around cylin- 
ders. We present solutions for Stokes flow past a station- 
ary cylinder (or, from another perspective, flow around 
a translating cylinder), the flow generated by a cylinder 
rotating at a constant angular velocity, and flow past a 
rotating cylinder (or, flow around a cylinder that is both 
rotating and translating). For each of these problems we 
have looked at the improvement in accuracy of the solu- 
tion as the finite element mesh is refined. 



4.1 Stokes Flow Past a Stationary Cylinder 

An exact solution exists for Stokes flow around a 
cylinder. The cylinder is surrounded by a fluid of in- 
finite extent. Far upstream of the cylinder the flow ap- 
proaches a constant one-dimensional velocity. The 
stream function for this flow is given by [12]: 
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Velocity components follow from the x and y derivatives 
of the stream function: 
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We created a regular, 1250-node (50 X 25) finite ele- 
ment mesh for this problem, and placed a particle at the 
center of the mesh. The flow was set equal to a constant 
value ([/„) at the boundaries of the mesh. The Reynolds 
number based on U„ was Re = 0.\. The resulting flow 
field is shown in Fig. 3. The same application was run 
with a 2450-node (70 X 35) mesh and a 3200-node 
(80 X 40) mesh. The root-mean-square (RMS) errors in 
the numerical solutions for the three meshes are, respec- 
tively, 5.059 %, 4.841 %, and 3.022 %. 

Separate errors were computed for each velocity 
component to better gauge the accuracy of the tech- 
nique: this gives information on the accuracy of the 
directions as well as magnitudes of the velocity vectors. 
As expected, the errors decrease as the mesh is refined. 
Errors for the x -component in the three meshes were 
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Fig. 3. Flow over a stationary cylinder — 1250-node mesh. 
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5.536 %, 5.127 %, and 3.070 %, and for the j-compo- 
nent, 2.514 %, 1.992 %, and 2.268 %. The 1250-node 
mesh is shown in Fig. 4, with the particle used in this 
simulation placed at the center. 

To test our moving-boundary method, we ran the 
same simulation in the frame of reference in which the 
particle is moving at a constant linear velocity, — C/„, 
and the fluid is stationary. The resulting flow is shown 
in Fig. 5. This is exactly (within machine accuracy) the 
solution that is obtained by taking the velocity field from 
the previous simulation and subtracting the velocity C/„. 
In this frame of reference it is easier to visualize how the 
particle moves the fluid out of its forward path, making 
it circulate above and below the cylinder, pulling fluid in 
toward itself at the rear. 

4.2 Stokes Flow Generated by a Rotating Cylinder 

To test the rotational component of our model, we 
analyzed the axisymmetric flow generated by a right. 



circular cylinder rotating at a constant angular velocity. 
In the ideal case, the fluid velocity goes to zero asymp- 
totically as the radial distance from the axis of symme- 
try goes to infinity. We used the same three finite ele- 
ment meshes as were used to test the model for flow 
around a stationary cylinder. The exact solution to this 
flow problem is given by Batchelor [13]: 



[/ = ■ 



R 



r 

R ' 



(33) 



in which R is the radius of the cylinder and r is the radial 
distance from the center of the cylinder. Figure 6 shows 
the resulting flow in the 50 X 25-node mesh. We calcu- 
lated the RMS error of the numerical solution computed 
on each of the three meshes. The error decreased from 
0.24 % to 0. 15 % from the 1250-node mesh to the 3200- 
node mesh. 



Fig. 4. 1250-node finite element mesh with a particle. 
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Fig. 5. Moving particles in an otherwise stationary fluid. 
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Fig. 6. Rotating cylinder in an otherwise stationary fluid. 



4.3 For Comparison: Stokes Flows Computed on a 
Fitted Mesh 

To compare our approach to approaches relying on 
remeshing, we computed solutions to both of the prob- 
lems presented above on a more "traditional" finite 
element mesh, constructed to fit the fluid domain, in- 
cluding the internal boundary. This mesh is shown in 
Fig. 7. It has a hole in the center of it to represent the 
particle, and the nodes near the cylinder have been 
moved to accommodate the hole. This mesh is similar to 
the (50 X 25)-node regular mesh, but with 1297 nodes. 
The extra nodes were added in the process of accommo- 
dating the particle. The error in the calculation per- 
formed with the fitted mesh was 0.68 %, which is con- 
siderably greater than the 0.24 % error observed for the 
(50 X 25)-node mesh using our new method. This addi- 



tional error is probably due to the irregular element 
shapes needed to form the hole in the mesh fitted. It is 
important to note that at least some of this discrepancy 
could be eliminated by optimizing the fitted mesh. 
However, this brings up another important point: in or- 
der to generate an accurate solution with a fitted mesh, 
given the necessary distortion (compared to a regular 
grid) to account for the hole, takes both effort and exper- 
tise. In contrast, generating the mesh for this problem 
using our technique is relatively trivial. 

4.4 Stokes Flow Past a Rotating Cylinder 

Figure 8 shows the flow past a rotating cylinder, the 
sum of the two flow fields already discussed, computed 
on the 1250-node mesh using our new technique. The 




Fig. 7. Finite element mesh representing a particle by a hole. 



86 



Volume 103, Number 1, January-February 1998 

Journal of Research of the National Institute of Standards and Technology 



' ^ > f p^^i ^ i ^ ^ » ^ * ^^ ' 







Fig. 8. Rotating cylinder in an otherwise uniform flow field. 



flow accelerates as it goes over the cylinder in the direc- 
tion of the rotation. This is the same flow field that is 
obtained by summing the velocity fields of the previous 
applications computed on the same mesh. This is ex- 
pected, since Stokes-flow solutions are linear, and con- 
sequently, superimposable. 

4.5 Multi-particle System 

We computed a Navier-Stokes solution of flow 
through a channel around a number of particles (effec- 
tively cylinders, since we are working in two dimen- 
sions) placed rather randomly across the channel. The 
Reynolds number is about 1 , based on a particle diame- 
ter. This illustrates the ability to place multiple objects 
in the flow. We selected this problem to allow a qualita- 
tive comparison with Martys and Garboczi [14]. They 



computed the flow field in two dimensions around 
roughly 100 particles using a straightforward finite dif- 
ference approach on a fine grid. Our solution here ap- 
pears qualitatively consistent with theirs. Figure 9 shows 
the flow through a channel around 8 distributed parti- 
cles. The obvious result is that the dominant flow fol- 
lows the path of least resistance, through the largest 
channels. 

Multi-particle systems do not require complicated 
meshing schemes with our technique. We can straight- 
forwardly model many particles without changing the 
background mesh, subject to the constraint of how fine 
the background mesh must be to capture the flow be- 
tween the particles, and the corresponding limits of 
computer memory. The particles can be located close to 
one another without distorting the mesh. However, they 
can approach only as close as allowed by the fineness of 




Fig. 9. Multi-particle flow. 
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the mesh, so that sufficient nodes remain between the 
particles to resolve the flow. On the coarse mesh used 
here, the particles are as close as allowable to maintain 
a qualitatively reasonable solution. Note that we do not 
currently allow the particles to contact each other. We 
have not yet accounted for particle-particle momentum 
transfer. This is not a major modification, and is dis- 
cussed further in the Discussion and Conclusions section. 



5. Examples — Falling Cylinders 

We sought experimental data with which to validate 
our 2-D model. While ample literature is available on 
falling spheres or other particle shapes, there appears to 
be relatively little on multiple falling cylinders. We con- 
ducted a series of crude experiments observing up to 
three cylindrical metal rods falling through glycerin. 
The rods are initially oriented horizontal. This is a stable 
configuration in terms of orientation: the drag on each 
cylinder is such that it discourages excursions away from 
the horizontal. (To stray from the horizontal, one end 
must fall faster. However, this increases the drag on that 
end, providing the necessary restoring force to maintain 
the cylinder in a horizontal orientation.) Multiple cylin- 
ders were parallel to each other. The rods had length-to- 
radius ratios of 10 to minimize end effects. We con- 
ducted a few experiments with rods whose aspect ratios 



were 20 and saw little or no difference compared to 
cylinders with aspect ratios of 10. The viscosity of glyc- 
erin is large, ^a = 11.4 g/(cm s) at 23 °C, but it behaves 
as a Newtonian fluid at this temperature. The density of 
glycerin at this temperature is 1 .25 g/cm^. The cylinders 
had diameters of 0.48 cm and lengths of about 5.5 cm. 
The cylinders were cold-rolled steel, with a density of 
about 7.74 g/cm^. Their settling speeds approached 16 
cm/s. The resulting Reynolds number is about 10 to 20 
based on the diameter of the cylinders. 

5.1 One Falling Cylinder 

We simulated a single particle (cylinder), initially at 
rest, settling under the influence of gravity. The simula- 
tion is 2-D. The boundary conditions are no-slip (zero 
velocity) on the top and bottom, and zero shear on the 
sides. This corresponds to a cylinder settling in a broad 
tank with a lid in contact with the liquid. In the simula- 
tion, the cylinder starts at a depth of about 1/3 of the 
total height of the tank, and falls about 1/3 of the height. 
Figure 10 shows the velocity vector- field for the flow 
around a single falling cylinder at a selected time step. 
The cylinder accelerates as it falls from its initial, sta- 
tionary position. The acceleration decreases during the 
simulation: the drag force on the cylinder increases as 
its speed increases, eventually balancing the force of 
gravity. 




Fig. 10. One falling cylinder. 
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An exact solution exists to determine the settling 
speed of a single cylinder falling in an incompressible 
fluid [15]. The normal and tangential stresses at the 
surface of the cylinder exert a drag force on the cylinder 
of magnitude D , given by the equation: 



in which C is a constant: 
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C-- 



In(^) 
where R is represented by: 
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The settling speed calculated from the exact solution is 
U= 16.85 cm/s. 

The simulation was run on meshes of several different 
grid spacings. The cylinder maintained a steady velocity 
only during the portion of the simulation during which 
it resided near the center of the mesh. Boundary effects 
were seen if the cylinder was too close to the upper or 
lower boundary. Our ability to simulate this situation 
was limited by the size of the mesh. The obvious cure 
for this problem is to translate the entire mesh at the 
settling speed of the particle, as was done by Johnson 
and Tezduyar [3]. On a domain with a grid spacing equal 
to 0.2 cm, on which each particle covers a 3 X 3-ele- 
ment area, the settling speed we computed was 15.9 
cm/s, which is in error by 5.6 % compared to the exact 
solution. On a mesh with a 0.12 cm grid-spacing, on 
which each particle covers a 5 X 5-element area, the 
settling speed was 17.0, which is in error by 0.89 %. 
The method appears to converge well. 

5.2 A Pair of Falling Cylinders 

An interesting result occurs experimentally when a 
pair of cylinders is placed several diameters apart near 
the top of the glycerol and released. The two cylinders 
fall at the same speed. They remain parallel and side-by- 
side, but the distance between them increases as they 
fall. The cylinders are observed to rotate about their 
horizontal, lengthwise axes in opposing directions. 
Looking at the cylinders end-on, the left one rotates 
clockwise, while the right one rotates counterclockwise. 
The explanation for this is as follows. The glycerol is 
quite viscous, so the momentum of each cylinder dif- 
fuses quickly toward the other. The two cylinders falling 
side-by-side restrict the flow between them, slowing the 
flow there and causing the pressure to increase between 



the cylinders. This causes some of the glycerin lying in 
the path between the cylinders to flow around the out- 
sides of the cylinders rather than between them, speed- 
ing the flow on the outsides of the cylinders. The pres- 
sure on the outsides of the cylinders is reduced, a 
manifestation of the Bernoulli effect. The resulting 
pressure gradient forces the cylinders apart. The rotation 
of the cylinders is the Magnus effect in reverse, where 
the pressure gradient (or outward diversion of the flow) 
causes the cylinders to counter-rotate. 

A finer mesh is required to reasonably simulate a pair 
of falling cylinders than is needed for the case of a single 
cylinder, because of the need to track the rotation of the 
particles with sufficient accuracy. Extraneous computa- 
tional constraints limited the size of the mesh we were 
able to use. The resolution of this problem is straightfor- 
ward, but is left for future work. Nonetheless, on the 
meshes that we used, the rotation we computed agrees 
qualitatively with what was observed experimentally. 
The cylinders fall side-by-side, moving apart and rotat- 
ing in opposite directions as they descend. Figure 11 
shows a snapshot of the velocity vector-field approxi- 
mately midway through the simulation of two cylinders 
falling. The horizontal acceleration is initially zero, and 
increases slowly at first. The horizontal velocity compo- 
nents in the figure are small and not easily visible. This 
is another case in which it would be helpful to translate 
the entire mesh at the speed of the center of mass of the 
cylinders, as was done by Johnson and Tezduyar [3]. 



6. Discussion and Conclusions 

We have created a novel technique for modeling two- 
phase flow of multiple particles moving in a fluid. An 
advantage of this approach is that the mesh does not 
distort as a result of the movement of the particles, 
obviating the need for remeshing and the associated 
numerical machinery. This technique is extensible to 
non-spherical particles. We have tested the technique by 
applying it to several simple linear problems for which 
exact solutions are known. We then illustrated the tech- 
nique for cylinders falling at low Reynolds numbers. 
Fluid velocities computed for locations very close to 
particles converge toward the exact solutions as the 
computational mesh is refined. 

An issue that has not yet been dealt with is the mod- 
eling of collisions between particles or with walls. For 
any technique this ultimately involves a judgment as to 
how close the particles may come before one considers 
them to contact. Physically, when the particles get too 
close, the assumptions of continuum mechanics are vio- 
lated. However, long before that, it is a practical neces- 
sity to limit how far the calculation proceeds, because as 
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Fig. 11. Two falling cylinders. 



the particles are allowed to approach more closely, the 
mesh must be correspondingly refined to capture the 
flow between the particles. The accuracy of our tech- 
nique is determined by the fineness of the background 
mesh. The close approach of particles creates a chal- 
lenge for remeshing schemes as well: frequent remesh- 
ing may be required to avoid unacceptable distortion as 
particles approach. We need to investigate the trade-off 
between the cost and accuracy of simulations as the 
background mesh is refined. If it turns out that too fine 
a mesh is required, a possibly useful compromise might 
be to use a hybrid scheme that uses a locally varying 
mesh, similar to what is used in other work, with the 
kind of particle definition we are using. In this way the 
mesh could be locally fine where necessary, but since 
our technique does not directly deform the mesh as the 
particles move it would limit the need for remeshing. We 
feel this technique is promising, and so have elected to 
present it at this time even though our numerical formu- 
lation needs to be refined to allow the finer and larger 
meshes required for more complex regimes of flow. 
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